package org.chorus;

import org.chorus.utils.KernelEstimator;

import jwo.landserf.structure.RasterMap;
import jwo.landserf.structure.SpatialModel;

//  **********************************************************************************
/** Represents the summary statistics describing a raster DEM.
 *  @author Jo Wood
 *  @version 1.3.2, 18th March, 2013.
 */
//  **********************************************************************************

/* This file is part of the Chorus river bed analysis software package. Chorus is free software:
 * you can redistribute it and/or modify it under the terms of the GNU General Public License
 * as published by the Free Software Foundation, either version 3 of the License, or (at your
 * option) any later version.
 * 
 * Chorus is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
 * without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  
 *  See the GNU General Public License for more details.
 *  
 * You should have received a copy of the GNU General Public License along with this
 * source code (see COPYING included with this source code). If not, see 
 * http://www.gnu.org/licenses/.
 */

public class RasterSummary 
{
	// ----------------------------- Object variables ------------------------------

	private RasterMap dem;
	private int numSamples;
	private double mean,stdev,skew,kurt;
	private float minVal,maxVal; 
	private KernelEstimator estimator;

	public RasterSummary(RasterMap dem)
	{
		this.dem = dem;		// Store underlying raster from which summary stats are found.
		calcSummaryStats();
	}

	// ---------------------------------- Methods ----------------------------------

	/** Calculates summary statistics of the DEM including mean, range, standard deviation,
	 *  skewness and kurtosis. 
	 */
	public void calcSummaryStats()
	{			
		// One pass calculation of mean, standard deviation and higher moments from Wikipedia
		// http://en.wikipedia.org/wiki/Algorithms_for_calculating_variance
		numSamples = 0;
		mean = 0;
		minVal = Float.MAX_VALUE;
		maxVal = -Float.MAX_VALUE;

		double M2 = 0;
		double M3 = 0;
		double M4 = 0;
		double delta,delta_n, delta_n2,term1;

		for (int row=0; row<dem.getNumRows();row++)
		{
			for (int col=0; col<dem.getNumCols();col++)
			{
				float elev = dem.getAttribute(row,col);

				if (elev != SpatialModel.NO_VALUE)
				{	
					minVal = Math.min(minVal,elev);
					maxVal = Math.max(maxVal, elev);
					numSamples++;
					delta = elev-mean;
					delta_n = delta / numSamples;
					delta_n2 = delta_n*delta_n;
					term1 = delta*delta_n*(numSamples-1);
					mean =  mean + delta_n;
					M4 = M4 + term1*delta_n2*(numSamples*numSamples - 3*numSamples + 3) + 6*delta_n2*M2 - 4*delta_n*M3;
					M3 = M3 + term1 * delta_n*(numSamples-2) - 3*delta_n*M2;
					M2 = M2 + term1;
				}
			}
		}

		stdev = Math.sqrt(M2/numSamples);
		skew = (Math.sqrt(numSamples)*M3)/Math.pow(M2, 1.5);
		kurt = (numSamples*M4)/(M2*M2) - 3.0;
	}
	
	/** Calculates variance of the given DEM in a single pass.
	 *  @param dem DEM to process.
	 */
	public static float calcVariance(RasterMap dem)
	{			
		// One pass calculation of mean, standard deviation and higher moments from Wikipedia
		// http://en.wikipedia.org/wiki/Algorithms_for_calculating_variance
		float numSamples = 0;
		double mean = 0;
		float minVal = Float.MAX_VALUE;
		float maxVal = -Float.MAX_VALUE;

		double M2 = 0;
		double delta,delta_n,term1;

		for (int row=0; row<dem.getNumRows();row++)
		{
			for (int col=0; col<dem.getNumCols();col++)
			{
				float elev = dem.getAttribute(row,col);

				if (elev != SpatialModel.NO_VALUE)
				{	
					minVal = Math.min(minVal,elev);
					maxVal = Math.max(maxVal, elev);
					numSamples++;
					delta = elev-mean;
					delta_n = delta / numSamples;
					term1 = delta*delta_n*(numSamples-1);
					mean =  mean + delta_n;
					M2 = M2 + term1;
				}
			}
		}
		return (float)(M2/numSamples);
	}

	/** Builds a Probability Density Function from the DEM values. Uses a Gaussian kernel with the
	 *  given bandwidth.
	 *  @param bandwidth Bandwidth of the KDE Gaussian function. Measured in the same units as the raster values.
	 */
	public void buildPDF(double bandwidth)
	{				
		estimator = new KernelEstimator(bandwidth,0.0000001);	

		// Select 50,000 random samples for PDF estimation.
		for (int i=0; i<50000; i++)
		{
			int row = (int)(Math.random()*dem.getNumRows());
			int col = (int)(Math.random()*dem.getNumCols());
			estimator.addValue(dem.getAttribute(row, col));
		}

		/* Sampling all cells appears to be too slow for large DEMs and unnecessary for PDF estimation.
		for (int row=0; row<dem.getNumRows();row++)
		{
			for (int col=0; col<dem.getNumCols();col++)
			{
				float elev = dem.getAttribute(row,col);

				if (elev != SpatialModel.NO_VALUE)
				{	
					estimator.addValue(elev);
				}
			}
		}
		 */

	}

	/** Reports the title of the raster being summarised.
	 *  @return Title of summarised raster.
	 */
	public String getTitle()
	{
		return dem.getHeader().getTitle();
	}

	/** Reports the mean of the raster values.
	 *  @return Mean of raster values.
	 */
	public float getMean()
	{
		return (float)mean;
	}

	/** Reports the standard deviation of the raster values.
	 *  @return Standard deviation of raster values.
	 */
	public float getStdev()
	{
		return (float)stdev;
	}

	/** Reports the skewness of the raster values.
	 *  @return Skewness of raster values.
	 */
	public float getSkew()
	{
		return (float)skew;
	}

	/** Reports the kurtosis of the raster values.
	 *  @return kurtosis of raster values.
	 */
	public float getKurt()
	{
		return (float)kurt;
	}

	/** Reports the minimum raster value.
	 *  @return Minimum of raster values.
	 */
	public float getMinVal()
	{
		return minVal;
	}

	/** Reports the maximum raster value.
	 *  @return Maximum of raster values.
	 */
	public float getMaxVal()
	{
		return maxVal;
	}

	/** Provides a probability estimate of the given value in the PDF.
	 *  @param x Value from which to estimate probability.
	 *  @return Estimated probability associated with the given x value of PDF.
	 */
	public float getProbability(double x)
	{
		return (float)(estimator.getProbability(x)/estimator.getPrecision());
	}


	/** Provides a textual representation of this summary.
	 *  @return Text representing this raster summary.
	 */
	public String toString()
	{
		StringBuffer text = new StringBuffer(dem.getHeader().getTitle()+"\n");
		text.append("Mean: "+mean+"\n");
		text.append("Stdev: "+stdev+"\n");
		text.append("Skewness: "+skew+"\n");
		text.append("Kurtosis: "+kurt+"\n");
		text.append("Min: "+minVal+"\n");
		text.append("Max: "+maxVal+"\n");
		return text.toString();
	}
}
